The samplr package provides a number of functions for sampling from continuous distributions.
Rejection sampling is a type of exact simulation method in numerical analysis and works for any distribution in \(\rm I\!R^m\) with a density. The idea of rejection sampling is to uniformly randomly sample the support and keep samples in the region under the graph of the density function.
Specifically, the steps of rejection sampling are:
Rejection sampling has two main assumptions:
The projectq2a function generates random deviates from a continuous random variable. The user supplies the probability density function of the continuous random variable, and the function utilizing rejection sampling to generate the deviates.
Here is an example of using projectq2a to sample from the standard uniform distribution.
The uniform distribution has probability density function:
\[\text{pdf = }\begin{cases}\frac{1}{b - a} \text{for } x \in [a, b]\\0 \text{ otherwise}\end{cases}\] We first generate a bunch of deviates and then plot their density in teal Then the true density is overlayed as a yellow line.
df_uniform_samples <-
data.frame(u = projectq2a(n = 10000,
pdf = dunif,
a = 0,
b = 1,
C = 1,
min = 0,
max = 1))
ggplot(df_uniform_samples, aes(x = u)) +
geom_density(color = NA, fill = "#35a09c") +
stat_function(aes(x = 0), fun = dunif, args = list(min = 0, max = 1), color = "#fff735", size = 1) +
labs(x = "X",
y = "Density") +
theme_minimal()Here is an example of using projectq2a to sample from a beta distribution with both shape parameters equal to 2.
The beta distribution has probability density function:
\[\text{pdf = }\begin{cases}\frac{x^{\alpha - 1}(1 - x)^{\beta - 1}}{B(\alpha, \beta)} \text{ for }0 \le x \le 1\\0 \text{ otherwise}\end{cases}\]
\[\text{where }B(\alpha, \beta) = \frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha + \beta)}\text{ with }\alpha, \beta \gt 0\]
We first generate a bunch of deviates and then plot their density in teal. Then the true density is overlayed as a yellow line.
df_beta_samples <-
data.frame(b = projectq2a(n = 10000,
pdf = dbeta,
a = 0,
b = 1,
C = 1.5,
shape1 = 2,
shape2 = 2))
ggplot(df_beta_samples, aes(x = b)) +
geom_density(color = NA, fill = "#35a09c") +
stat_function(aes(x = 0), fun = dbeta, args = list(shape1 = 2, shape2 = 2), color = "#fff735", size = 1) +
labs(x = "X",
y = "Density") +
theme_minimal()Here is an example of using projectq2a to sample from a custom distribution.
In this case, the custom density function is that of a beta with shape parameters 2 and 5.
\[\text{pdf = }\begin{cases}30x(1 - x)^4 \text{ for }0 \le x \le 1\\0 \text{ otherwise}\end{cases}\]
dcustom <- function (x) {
sapply(X = x,
FUN = function(x) {
ifelse(0 < x & x < 1,
(gamma(2+5)/(gamma(2)*gamma(5)))*(x^(2-1)*(1-x)^(5-1)),
0)
})
}We first generate a bunch of deviates and then plot their density in teal. Then the true density is overlayed as a yellow line.
df_custom_samples <-
data.frame(b = projectq2a(n = 10000,
pdf = dcustom,
a = 0,
b = 1,
C = 2.5))
ggplot(df_custom_samples, aes(x = b)) +
geom_density(color = NA, fill = "#35a09c") +
stat_function(aes(x = 0), fun = dcustom, color = "#fff735", size = 1) +
labs(x = "X",
y = "Density") +
theme_minimal()The projectq3a function generates random deviates from a continuous 2D distribution defined on a square. The user supplies the probability density function, and the function utilizes rejection sampling to generate the deviate pairs.
We first demo sampling from a uniform distribution defined on the unit square. In order to do so, we must first define a joint probability density function which takes x and y as arguments for the two variables.
\[\text{pdf = }\begin{cases}\frac{1}{(b - a)^2} \text{for } x, y \in [a, b]\\0 \text{ otherwise}\end{cases}\]
d2dunif <- function(x, y, min = 0, max = 1) {
if(min <= x && x <= max && min <= y && y <= max)
(max - min)^(-2)
else
0
}Then we could examine what the distribution looks like in 3D by computing the probability densities for each x and y pair in a systematic sampling across the support.
l <- list(x = seq(0, 1, 0.5),
y = seq(0, 1, 0.5))
z <- matrix(rep(NA, length(l$x) * length(l$y)),
nrow = length(l$x))
for(r in 1:nrow(z)) {
for(c in 1:ncol(z)) {
z[r, c] <- d2dunif(x = l$x[r], y = l$y[c])
}
}
l$z <- z
plot_ly(x = l$x, y = l$y, z = l$z) %>%
add_surface() %>%
layout(scene = list(xaxis = list(title = "X"),
yaxis = list(title = "Y"),
zaxis = list(title = "Density")))Finally, we can sample from the joint distribution utilizing rejection sampling and estimate the kernel density.
We first demo sampling from a beta distribution defined on the unit square. In order to do so, we must first define a joint probability density function which takes x and y as arguments for the two variables.
\[\text{pdf = }\begin{cases}\frac{x^{\alpha - 1}(1 - x)^{\beta - 1} + y^{\alpha - 1}(1 - y)^{\beta - 1}}{2B(\alpha, \beta)}\text{ for }0 \le x, y \le 1\\0\text{ otherwise}\end{cases}\]
\[\text{where }B(\alpha, \beta) = \frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha + \beta)}, \text{ with }\alpha, \beta \gt 0\]
d2dbeta <- function(x, y, shape1, shape2) {
if(0 <= x && x <= 1 && 0 <= y && y <= 1)
(dbeta(x = x, shape1 = shape1, shape2 = shape2) +
dbeta(x = y, shape1 = shape1, shape2 = shape2)) / 2
else
0
}Then we could examine what the distribution looks like in 3D by computing the probability densities for each x and y pair in a systematic sampling across the support.
l <- list(x = seq(0, 1, 0.01),
y = seq(0, 1, 0.01))
z <- matrix(rep(NA, length(l$x) * length(l$y)),
nrow = length(l$x))
for(r in 1:nrow(z)) {
for(c in 1:ncol(z)) {
z[r, c] <- d2dbeta(x = l$x[r], y = l$y[c], shape1 = 2, shape2 = 2)
}
}
l$z <- z
plot_ly(x = l$x, y = l$y, z = l$z) %>%
add_surface() %>%
layout(scene = list(xaxis = list(title = "X"),
yaxis = list(title = "Y"),
zaxis = list(title = "Density")))Finally, we can sample from the joint distribution utilizing rejection sampling and estimate the kernel density.
We next demo sampling from a custom joint distribution defined on the square from -1 to +1. Again, we first define the joint pdf function which takes x and y as quantile arguments.
\[\text{pdf = }\begin{cases}\frac{3}{8}(x^2 + y^2) \text{for } x, y \in [-1, +1]\\0 \text{ otherwise}\end{cases}\]
d2dcirclecontour <- function(x, y) {
if(-1 <= x && x <= 1 && -1 <= y && y <= 1)
(3/8)*(x^2 + y^2)
else
0
}Then we could examine what the distribution looks like in 3D by computing the probability densities for each x and y pair in a systematic sampling across the support.
l <- list(x = seq(-1, 1, 0.1),
y = seq(-1, 1, 0.1))
z <- matrix(rep(NA, length(l$x) * length(l$y)),
nrow = length(l$x))
for(r in 1:nrow(z)) {
for(c in 1:ncol(z)) {
z[r, c] <- d2dcirclecontour(x = l$x[r], y = l$y[c])
}
}
l$z <- z
plot_ly(x = l$x, y = l$y, z = l$z) %>%
add_surface() %>%
layout(scene = list(xaxis = list(title = "X"),
yaxis = list(title = "Y"),
zaxis = list(title = "Density")))Finally, we can sample from the joint distribution utilizing rejection sampling and estimate the kernel density.